clear all
close all
clc

K = 30; % K is the length of decoded bits
S = K+24;             % Message length in bits, including CRC, K > 30
L = 8; % List length, a power of two, [1 2 4 8]
E = 128; % polar码输出序列长度
M = 8; % 天线数量
K_R = 3; % Rician-K-factor

crcLen = 24;        % Number of CRC bits for DL, Section 5.1, [6]
% [DFT,angle] = Function_ULA_code_DFT(M); % 设置DFT矩阵
x = 0 : (M - 1);
SNR_dB = 10:5:30; 
SNR = 10.^(SNR_dB/10);
g_sym = 0.43 : 0.01 : 0.51;
ii_max = 500;

theta_b = 100; % 用户位置
theta_e_s = 93 ; % 窃听者位置

[DFT,angle] = Function_ULA_code_DFT2(M,theta_b); % 设置DFT矩阵

% angleTemp(id_SNR,ii_gs,ii_it) = zeros(5,9,ii_max);
% g_sym_fin(id_SNR,ii_it) = zeros(5,ii_max);
% angle_d_fin(id_SNR,ii_it) = zeros(5,ii_max);
% numBitErrs_B(id_SNR,ii_it) = zeros(5,ii_max);
% numBitErrs_E(id_SNR,ii_it) = zeros(5,ii_max);

comb = zeros(73,1);

trblk = randi([0 1],K,1,'int8'); % 生成初始序列
dmrsOri = randi([0,1],[18,1]); % 生成导频序列
an_i = randi([0 1],146,1,'int8'); % 生成AN序列
    
msgcrc = nrCRCEncode(cast(trblk,'int8'),'24C');
encOut = nrPolarEncode(msgcrc,E); % Polar码
    
sym = nrSymbolModulate(encOut,'QPSK'); % QPSK调制
dmrs = nrSymbolModulate(dmrsOri,'QPSK'); % QPSK调制
an_q = nrSymbolModulate(an_i,'QPSK'); % QPSK调制

t = 1;
j = 1;
for ii = 1:73
    if mod(ii,8) > 0
        comb(ii,1) =  sym(t);
        t = t + 1;
    else
        comb(ii,1) =  dmrs(j);
        j = j + 1;
    end
end
P_T = mean(abs(comb).^2);
    
for thetaE_indice = 1 : 1
    ini_b = 180;
    ini_e = 180;
    for ii = 1:M
       delta_b = abs(theta_b - angle(ii));
       if delta_b < ini_b
            ini_b = delta_b;
            DFTindice = ii;
       end
    end
        
    for ii = 1:M
       delta_e = abs(theta_e_s(thetaE_indice) - angle(ii));
        if delta_e < ini_e && ii ~= DFTindice
            ini_e = delta_e;
            ANindice = ii;
        end
    end
        
    g = 0 :(M-1);
    sym_beam = exp(-1j * pi * g * cos(angle(DFTindice) / 180 * pi))/sqrt(M);
%     sym_beam2 = (DFT(:,DFTindice)/sqrt(M)).';   

    for id_SNR = 1:5
        SINR_e_max = 0;
        % 生成合法用户的莱斯信道
        H_0b = exp(1j * pi * x * cos(theta_b / 180 * pi)); % 生成莱斯信道的LOS成分
        H_rb = randn(1,M); % 生成莱斯信道的NLos成分
        H_b = (sqrt(K_R / (1 + K_R)) * H_0b + sqrt(1 / (1 + K_R)) * H_rb); % 生成莱斯信道

        % 生成窃听者的莱斯信道
        H_0e = exp(1j * pi * x * cos(theta_e_s(thetaE_indice) / 180 * pi)); % 生成莱斯信道的LOS成分
        H_re = randn(1,M); % 生成莱斯信道的NLos成分
        H_e = (sqrt(K_R / (1 + K_R)) * H_0e + sqrt(1 / (1 + K_R)) * H_re); % 生成莱斯信道
            
        noiseVar(id_SNR) = 1/SNR(id_SNR);
        noise_B = sqrt(noiseVar(id_SNR)/2)*(randn(73,1)+1i*randn(73,1));
        noise_E = sqrt(noiseVar(id_SNR)/2)*(randn(73,1)+1i*randn(73,1));
        for ii_gs = 1 : 9  
            SINR = 0;
            mark = 0;
            angleTemp(id_SNR,ii_gs) = angle(ANindice) - 0.1;
            while SINR < 0.001 % P(sinr_e < 40) <= 0.02
                angleTemp(id_SNR,ii_gs) = angleTemp(id_SNR,ii_gs) + 0.1;
                an_beam = exp(-1j * pi * g * cos(angleTemp(id_SNR,ii_gs) / 180 * pi))/sqrt(M);
                    
                % 计算用户SINR的概率密度函数
                K_slide_b(id_SNR) = abs(sym_beam * H_0b.').^2 * K_R;
                K_slide_an(id_SNR) = abs(an_beam * H_0b.').^2 * K_R;
                m_b_slide(id_SNR) = (K_slide_b(id_SNR)+1).^2 / (2*K_slide_b(id_SNR) + 1);
                m_an_slide(id_SNR) = (K_slide_an(id_SNR)+1).^2 / (2*K_slide_an(id_SNR) + 1);
                r_slide_b(id_SNR) = (g_sym(ii_gs).^2) / noiseVar(id_SNR);
                r_slide_an(id_SNR) = ((1 - g_sym(ii_gs)).^2) / noiseVar(id_SNR);
                r_slide_b_ba(id_SNR) = ((K_slide_b(id_SNR)+1) * r_slide_b(id_SNR)) / (K_R + 1);
                r_slide_an_ba(id_SNR) = ((K_slide_an(id_SNR)+1) * r_slide_an(id_SNR)) / (K_R + 1);
                C_Ab = (m_b_slide(id_SNR) ./ r_slide_b_ba(id_SNR)).^m_b_slide(id_SNR);
                C_An = (m_an_slide(id_SNR) ./ r_slide_an_ba(id_SNR)).^m_an_slide(id_SNR);
                C_Gb = gamma(m_b_slide(id_SNR));
                C_Gn = gamma(m_an_slide(id_SNR)); % gamma(m2)
                C_eb = (-1) * m_b_slide(id_SNR) ./ r_slide_b_ba(id_SNR);
                C_en = (-1) * m_an_slide(id_SNR) ./ r_slide_an_ba(id_SNR);
                m1 = m_b_slide(id_SNR);
                m2 = m_an_slide(id_SNR);
                C_all = (C_Ab .* C_An) / (C_Gb .* C_Gn) .* (1 / noiseVar(id_SNR)).^(m1 + m2).* exp((-1) * C_en) .* (noiseVar(id_SNR).^(m1 + m2));
                b = m1 + m2 + 1;
                f_sinr = @(y) C_all .* gamma(m2) .* kummerU(m2, b, ((-1) .* (C_eb .* y + C_en))) .* y.^(m1 - 1) .* exp((-1) .* ((-1) .* (C_eb .* y + C_en)));
                SINR = quadgk(f_sinr,0,1.4);        
                    
                if SINR < 0.003
                   mark = 1; % 至少有一组解满足SINR_b的阈值条件
                end
            end
            if mark == 1
                an_temp = angleTemp(id_SNR,ii_gs) - 0.2;
                an_beam = exp(-1j * pi * g * cos(an_temp / 180 * pi))/sqrt(M);
            end
%             trans = g_sym(ii_gs) * comb * sym_beam + (1 - g_sym(ii_gs)) * an_q * an_beam; % 使用DFT矩阵的某一列作为波束
%             receiveB1 = trans * H_b.';
%             receiveE1 = trans * H_e.';
%                 
%             receiveB2 = receiveB1 + noise_B;
%             receiveE2 = receiveE1 + noise_E;
% 
%             for ii = 1 : 9
%                 HP_E(ii) = receiveE2(8 * ii) ./ comb(8 * ii);
%             end
%                 
%             HR_E = interp1(1:8:64,HP_E(1:8),1:1:64,'pchip');
% 
%             t = 1;
%             for ii = 1:73
%                if mod(ii,8) > 0
%                    Rsym_B(t,1) =  receiveB2(ii);
%                    Rsym_E(t,1) =  receiveE2(ii);
%                    t = t + 1;
%                end
%             end
% 
%             demod_B = nrSymbolDemodulate(Rsym_B,'QPSK',noiseVar(id_SNR)); % QPSK解调
%             demod_E = nrSymbolDemodulate(Rsym_E ./ HR_E.','QPSK',noiseVar(id_SNR)); % QPSK解调
%                 
%             decBlk_B2 = nrPolarDecode(demod_B,S,E,L); % Polar解码
%             decBlk_E2 = nrPolarDecode(demod_E,S,E,L); % Polar解码
%                 
%             [decBlk_B,errFlag_B] = nrCRCDecode(decBlk_B2,'24C');
%             [decBlk_E,errFlag_E] = nrCRCDecode(decBlk_E2,'24C');
%                 
%             % 计算误码率，并进行比较
%             numBitErrs_E_temp = biterr(decBlk_E,trblk) / K;

            % 计算窃听者SINR的概率密度函数
            K_slide_b(id_SNR) = abs(sym_beam * H_0e.').^2 * K_R;
            K_slide_an(id_SNR) = abs(an_beam * H_0e.').^2 * K_R;
            m_b_slide(id_SNR) = (K_slide_b(id_SNR)+1).^2 / (2*K_slide_b(id_SNR) + 1);
            m_an_slide(id_SNR) = (K_slide_an(id_SNR)+1).^2 / (2*K_slide_an(id_SNR) + 1);
            r_slide_b(id_SNR) = (g_sym(ii_gs).^2) / noiseVar(id_SNR);
            r_slide_an(id_SNR) = ((1 - g_sym(ii_gs)).^2) / noiseVar(id_SNR);
            r_slide_b_ba(id_SNR) = ((K_slide_b(id_SNR)+1) * r_slide_b(id_SNR)) / (K_R + 1);
            r_slide_an_ba(id_SNR) = ((K_slide_an(id_SNR)+1) * r_slide_an(id_SNR)) / (K_R + 1);

            C_Ab = (m_b_slide(id_SNR) ./ r_slide_b_ba(id_SNR)).^m_b_slide(id_SNR);
            C_An = (m_an_slide(id_SNR) ./ r_slide_an_ba(id_SNR)).^m_an_slide(id_SNR);
            C_Gb = gamma(m_b_slide(id_SNR));
            C_Gn = gamma(m_an_slide(id_SNR)); % gamma(m2)
            C_eb = (-1) * m_b_slide(id_SNR) ./ r_slide_b_ba(id_SNR);
            C_en = (-1) * m_an_slide(id_SNR) ./ r_slide_an_ba(id_SNR);
            m1 = m_b_slide(id_SNR);
            m2 = m_an_slide(id_SNR);
            C_all = (C_Ab .* C_An) / (C_Gb .* C_Gn) .* (1 / noiseVar(id_SNR)).^(m1 + m2).* exp((-1) * C_en) .* (noiseVar(id_SNR).^(m1 + m2));
            b = m1 + m2 + 1;
            f_sinr = @(y) C_all .* gamma(m2) .* kummerU(m2, b, ((-1) .* (C_eb .* y + C_en))) .* y.^(m1 - 1) .* exp((-1) .* ((-1) .* (C_eb .* y + C_en)));
                   
            SINR_e_temp = quadgk(f_sinr,0,2);

            if SINR_e_temp > SINR_e_max && mark == 1
                g_sym_fin(id_SNR) = g_sym(ii_gs);
                angle_d_fin(id_SNR) = an_temp;
                SINR_e_max = SINR_e_temp;
            end   
        end
    end
end

% 运行
for ii_it = 1:ii_max
    
    trblk = randi([0 1],K,1,'int8'); % 生成初始序列
    dmrsOri = randi([0,1],[18,1]); % 生成导频序列
    an_i = randi([0 1],146,1,'int8'); % 生成AN序列
    
    msgcrc = nrCRCEncode(cast(trblk,'int8'),'24C');
    encOut = nrPolarEncode(msgcrc,E); % Polar码
    
    sym = nrSymbolModulate(encOut,'QPSK'); % QPSK调制
    dmrs = nrSymbolModulate(dmrsOri,'QPSK'); % QPSK调制
    an_q = nrSymbolModulate(an_i,'QPSK'); % QPSK调制

    t = 1;
    j = 1;
    for ii = 1:73
        if mod(ii,8) > 0
            comb(ii,1) =  sym(t);
            t = t + 1;
        else
            comb(ii,1) =  dmrs(j);
            j = j + 1;
        end
    end
    P_T = mean(abs(comb).^2);
    
    for thetaE_indice = 1 : 1
        ini_b = 180;
        ini_e = 180;
        for ii = 1:M
           delta_b = abs(theta_b - angle(ii));
           if delta_b < ini_b
                ini_b = delta_b;
                DFTindice = ii;
           end
        end
        
        g = 0 :(M-1);
        sym_beam = exp(-1j * pi * g * cos(angle(DFTindice) / 180 * pi))/sqrt(M);
    
        for id_SNR = 1:5    
            % 生成合法用户的莱斯信道
            H_0b = exp(1j * pi * x * cos(theta_b / 180 * pi)); % 生成莱斯信道的LOS成分
            H_rb = randn(1,M); % 生成莱斯信道的NLos成分
            H_b = (sqrt(K_R / (1 + K_R)) * H_0b + sqrt(1 / (1 + K_R)) * H_rb); % 生成莱斯信道

            % 生成窃听者的莱斯信道
            H_0e = exp(1j * pi * x * cos(theta_e_s(thetaE_indice) / 180 * pi)); % 生成莱斯信道的LOS成分
            H_re = randn(1,M); % 生成莱斯信道的NLos成分
            H_e = (sqrt(K_R / (1 + K_R)) * H_0e + sqrt(1 / (1 + K_R)) * H_re); % 生成莱斯信道
            
            noiseVar(id_SNR) = 1/SNR(id_SNR);
            noise_B = sqrt(noiseVar(id_SNR)/2)*(randn(73,1)+1i*randn(73,1));
            noise_E = sqrt(noiseVar(id_SNR)/2)*(randn(73,1)+1i*randn(73,1));

            an_beam = exp(-1j * pi * g * cos(angle_d_fin(id_SNR) / 180 * pi))/sqrt(M);
            g_sym = g_sym_fin(id_SNR);
            
            trans = sqrt(g_sym) * comb * sym_beam + sqrt(1 - g_sym) * an_q * an_beam; % 使用DFT矩阵的某一列作为波束
            receiveB1 = trans * H_b.';
            receiveE1 = trans * H_e.';
                
            receiveB2 = receiveB1 + noise_B;
            receiveE2 = receiveE1 + noise_E;
            
            for ii = 1 : 9
                HP_E(ii) = receiveE2(8 * ii) ./ comb(8 * ii);
            end
                
            HR_E = interp1(1:8:64,HP_E(1:8),1:1:64,'pchip');

            t = 1;
            for ii = 1:73
               if mod(ii,8) > 0
                   Rsym_B(t,1) =  receiveB2(ii);
                   Rsym_E(t,1) =  receiveE2(ii);
                   t = t + 1;
               end
            end

            demod_B = nrSymbolDemodulate(Rsym_B,'QPSK',noiseVar(id_SNR)); % QPSK解调
            demod_E = nrSymbolDemodulate(Rsym_E ./ HR_E.','QPSK',noiseVar(id_SNR)); % QPSK解调
            
            demod_B_test = nrSymbolDemodulate(Rsym_B,'QPSK','DecisionType','hard'); % QPSK解调
            demod_E_test = nrSymbolDemodulate(Rsym_E ./ HR_E.','QPSK','DecisionType','hard'); % QPSK解调
            
            SER_B(id_SNR,ii_it) = biterr(demod_B_test,encOut) / E;     
            SER_E(id_SNR,ii_it) = biterr(demod_E_test,encOut) / E;  
            
            decBlk_B2 = nrPolarDecode(demod_B,S,E,L); % Polar解码
            decBlk_E2 = nrPolarDecode(demod_E,S,E,L); % Polar解码
                
            [decBlk_B,errFlag_B] = nrCRCDecode(decBlk_B2,'24C');
            [decBlk_E,errFlag_E] = nrCRCDecode(decBlk_E2,'24C');            
                          
            numBitErrs_B(id_SNR,ii_it) = biterr(decBlk_B,trblk) / K;
            numBitErrs_E(id_SNR,ii_it) = biterr(decBlk_E,trblk) / K;
            if numBitErrs_B(id_SNR,ii_it) ~= 0
                numBlockErrs_B(id_SNR,ii_it) = 1;
            else
                numBlockErrs_B(id_SNR,ii_it) = 0;
            end
            if numBitErrs_E(id_SNR,ii_it) ~= 0
                numBlockErrs_E(id_SNR,ii_it) = 1;
            else
                numBlockErrs_E(id_SNR,ii_it) = 0;
            end   
        end
    end       
end

b1 = mean(numBlockErrs_B,2);
b2 = mean(numBlockErrs_E,2);

s1 = mean(SER_B,2);
s2 = mean(SER_E,2);

figure,
plot(SNR_dB,s1,'*-','MarkerEdgeColor','R','linewidth',1);
xlabel('SNR-dB'),ylabel('SER')
hold on
plot(SNR_dB,s2,'*-','MarkerEdgeColor','B','linewidth',1);
title('添加AN但有角度偏移时的SER性能 (Bob at 100, Eve at 93)')
legend('Bob','Eve')

figure,
plot(SNR_dB,b1,'*-','MarkerEdgeColor','R','linewidth',1);
xlabel('SNR-dB'),ylabel('BLER')
hold on
plot(SNR_dB,b2,'*-','MarkerEdgeColor','B','linewidth',1);
title('添加AN但有角度偏移时的BLER性能 (Bob at 100, Eve at 93)')
legend('Bob','Eve')










